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In the Cauchy problem of general relativity one considers initial data that satisfies certain con- 
straints. The evolution equations guarantee that the evolved variables will satisfy the constraints 
at later instants of time. This is only true within the domain of dependence of the initial data. 
^SJ ' If one wishes to consider situations where the evolutions are studied for longer intervals than the 

, size of the domain of dependence, as is usually the case in three dimensional numerical relativity, 

■ one needs to give boundary data. The boundary data should be specified in such a way that the 

constraints are satisfied everywhere, at all times. In this paper we address this problem for the case 
of general relativity linearized around Minkowski space using the generalized Einstein-Christoffel 
^ . symmetric hyperbolic system of evolution equations. We study the evolution equations for the con- 

' straints, specify boundary conditions for them that make them well posed and further choose these 

boundary conditions in such a way that the evolution equations for the metric variables are also 
well posed. We also consider the case of a manifold with a non-smooth boundary, as is the usual 
case of the cubic boxes commonly used in numerical relativity. The techniques discussed should be 
applicable to more general cases, as linearizations around more complicated backgrounds, and may 
be used to establish well posedness in the full non-linear case. 
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I. INTRODUCTION 



o 
a^ 
o 

In the Cauchy problem of general relativity, Einstein's equations split into a set of evolution equations and a set of 
constraint equations. Given initial data that satisfies the constraints one can show that the evolution equations imply 
O ■ that the constraints continue to hold in the domain of dependence of the initial slice. However, in numerical relativity, 
'■ due to finite computer resources, one usually solves the field equations in a domain of the form {t, x^) £ TZ x with 
$H . ^ sl bounded space-like surface with boundary dil. In this case the domain of dependence of the initial slice is "too 
small" and one wants to evolve beyond it. Then the question that arises is what boundary conditions to give at the 
artificial boundary TZ x dCt . This problem has two aspects to it. First, one wishes to prescribe boundary conditions 
, . ■ that make the problem a well posed one which preserves the constraints. In addition to that, one might desire to 
rS I embed in the boundary conditions some physically appealing property (for instance that no gravitational radiation 
. enters the domain). These two aspects are in principle separate. In this paper we will concentrate on the first one, 
namely, how to prescribe consistent boundary conditions. The construction we propose ends up giving "Dirichlet 
or Neumann-like" boundary conditions on some components of the metric, with some free sources. The latter allow 
freedom to consider boundary conditions for any spacetime in any slicing. 

If the field equations are cast into a first order symmetric hyperbolic form there is a known prescription for 
achieving well posedness for an initial-boundary value problem (this is discussed in detail in section II), namely 
maximal dissipative boundary conditions Well posedness is a necessary condition for implementing a stable 
numerical code in the sense of Lax's theorem [^j. However, additional work is needed if one wishes to have a code 
that is not only stable but that preserves the constraints throughout the domain of evolution. As we stated above, 
this implies giving initial and boundary data that guarantees that the constraints are satisfied. Traditionally, most 
numerical relativity treatments have been careful to impose initial data that satisfies the constraints. However, very 
rarely boundary conditions that lead to well posedness are used, and much less frequently are they consistent with the 
constraints. Only recently, following Friedrich and Nagy has work with numerical relativity in mind started along 
these lines , Q , 1^ , , Q . Of particular interest is the recent paper of Szilagyi and Winicour , which outlines a 
construction with several points in common with the one we describe here. 

In this article we derive well posed, constraint-preserving boundary conditions for the generalized pO[ | Einstein- 
Christoffel (EC) [0 symmetric hyperbolic formulations of Einstein's equations, when linearized around flat spacetime. 
The procedure consists in studying the evolution system for the constraints and making sure that the corresponding 
initial-boundary value problem is well-posed through appropriate boundary conditions. There is some freedom in the 
specification of the latter. Next, we translate them into boundary conditions for the variables of the main evolution 
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system. The main difhculty consists in making use of the freedom that one has in the boundary conditions for the 
evolution system of the constraints in such a way that the resulting boundary conditions for the main evolution system 
ensure well-posedness. We consider the case of a non-smooth cubic boundary as is of interest in numerical relativity. 
This involves the additional complication of being careful about ensuring compatibility at the edges joining the faces. 
To our knowledge, this is the first detailed analysis of the initial-boundary value problem with a non-smooth boundary 
in the gravitational context. 

The organization of this paper is as follows: in section |ll| we review the basic technique of energy estimates used 
to prove well posedness. In section [II we review the generalized EC symmetric hyperbolic formulation, write down 
the evolution system for the constraints and analyze under what conditions the latter is symmetric hyperbolic. In 
section EV we compute the characteristic variables for the main evolution system and for the system that evolves the 
constraints. These characteristic variables are a key element for prescribing the boundary conditions that lead to 
well-posedness. In section ^ we write down the constraint preserving boundary conditions in terms of the variables 



of the main evolution system. In section VI we derive the necessary en ergy estimates to show that all the systems 



involved in our constraint-preserving treatment are well posed. In section VII we summarize the constraint-preserving 
construction of this paper. We end with a discussion of possible future improvements to and applications of the 
boundary conditions introduced in this paper. 



II. BASIC ENERGY ESTIMATES 



In this section we review some basic notions of energy estimates for systems of partial differential equations, as 
discussed, for instance, in Consider a first order in time and space linear evolution system of the form 

dtu = A^djU, (1) 

where u = u{t,x^) is a vector valued function, t > 0, S 51, and where the matrices , and A^ are constant. 
The symbol P{n) is defined as P := A^Ui. The system is symmetric (or symmetrizable) hyperbolic if there is a 
symmetrizer for P. That is, a positive definite, Hermitian matrix H independent of such that HP = P^H for all 
with Ej=i...3"? = 1- 

In order to show well-posedness of the initial value problem for such a system, one usually derives a bound for the 
energy 

E{t) = I {u,Hu)d^x , (2) 
Jn 

where (.,.) denotes the standard scalar product. Taking a time derivative of (||), using equation (|l|), the fact that 
HA^ are symmetric matrices, and Gauss' theorem, one obtains 

^E{t)^ [ 2{u,HA^dju)d^x = I dj{u,HA^u)d^x ^ [ {u,HP{n)u) da, 
dt Jn Jn Jan 

where n is now the unit outward normal to the boundary dil of the domain. For simplicity let's assume that P(n) has 
only the eigenvalues ±1 and 0. Let u^~\ and u^^^ denote the projections of u onto the eigenspaces corresponding 
to the eigenvalues 1,-1 and 0, respectively. That is, u — -t- w^^-* + u'^"^ and P{n)u — u^^'^ — u'^~\ Then 

{u,HP{n)u) = {u''+\Hu'-+'^) - {u^-\Hu''-^). 

If we impose boundary conditions of the form = Ru^^'' with R "small enough" so that R^HR < H it 
follows that E{t) < E{0) for all t > 0. These kind of energy estimates are a key ingredient in well-posedness proofs. 
One can generalize this result to boundary conditions of the form 

= +5, (3) 

where g = g{t, x^) is a prescribed function at the boundary and R^HR < H, as before. In this case, an energy 
estimate can be obtained as follows: We first choose a function -0 that satisfies V'^^^ = R'4^^~^ + g &t the boundary. 
Then, we consider the variable u = u — instead of u which now satisfies the homogeneous boundary condition 
= Ru^^^ and the modified evolution equation 



dtu = A^dju + F, 
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where F = A^djip — dtip is a forcing term. One now repeats the estimate for the energy defined in (^, with u replaced 
by u, and obtains 

-i-E{t) < 2 / {u,HF)d^x < 2\\u\\ ■ \\F\\, 

where ||w|| = E^^"^, \\F\\ = {J {F, HF)d^xY^^ and where we have used Schwarz's inequahty. Therefore, we obtain the 
estimate 

\\nit,.)\\<\m-)\\+ f\\Fis,.)\\ds. 

Jo 

Similar estimates can be obtained for systems with non-constant coefficients, as is the case, for instance, of lin- 
earizations around a given non- Minkowski metric, as in black hole perturbations. In the non-linear case the proofs of 
bounds are only for finite amounts of time, since solutions can blow up in a finite time starting from initial data with 
finite energy. 

The existence of an energy estimate implies that the initial-boundary value problem is well posed. By this we mean 
that there exists a unique smooth solution to the problem for which the energy estimate holds It should 

be noted that these results are in principle only valid for smooth boundaries. Later in this paper we will discuss 
boundaries that are non-smooth. It should be understood that in those cases the energy estimates we derive do not 
necessarily imply the existence and uniqueness of smooth solutions. 



III. THE FIELD EQUATIONS AND THEIR CONSTRAINT PROPAGATION 

In this section we present the field equations that we will use through this paper, and analyze the constraints 
propagation. In the first subsection we present the evolution equations for the main variables within the generalized 
Einstcin-Christoffel symmetric hyperbolic formulation of Einstein's equations. In the second subsection we analyze, 
in the fully non-linear case, the evolution equations for the constraints within this formulation and derive necessary 
conditions for the latter to be symmetrizable. We will need this system to be symmetric hyperbolic in order to 
derive an energy estimate in the way we sketched in Section |l[ Rather surprisingly, it turns out that in the original 
EC system the constraints' propagation does not seem to be symmetrizable. Imposing the symmetric hyperbolicity 
condition naturally restricts the free parameter of the system to an open interval. As a side note, there seems to 
be some correlation between the stability properties of the system found in numerical experiments and this natural 
choice. 



A. The formulation 



of evolution equations for the three-metric ), the extrinsic 
curvature {Kij ) , and some extra variables fkij that are introduced in order to make the system first order in space is 
derived: 

dQ9tj = -'2.K,j, (4) 

doK,, - -d'^fkrj+l.O., (5) 

dofki] = -dkKij + l.o. , (6) 

where d'^ = g'^^di and where Oq = {dt — £p)/N is the derivative operator along the normal to the spatial t = const. 
slices. The shift is an a priori prescribed function on spacetime while the lapse N is determined by TV = g^^^ e*^, 
with Q a priori prescribed. Here and in the following l.o. stands for "lower order terms". These terms depend on 
gij , Kij , fkij , Q and /3* but not on the derivatives of , Kij or f-^ij . Provided that the constraints (see equation (|^) 
below) are satisfied, the spatial derivatives, dfey = dugij, of the three-metric are obtained from 

du^J = 2 /fc,, + rj (/^.), ^ - f)^^}j + ^ 5y - /IJ , (7) 

where 7y is a free parameter with the only restriction rj ^ 0. 

The evolution system (^,||,||), besides being symmetric hyperbolic, has the additional feature that all characteristic 
speeds (with respect to the normal derivative operator do) are either or ±1, so that all characteristic modes lie 
either along the light cone or along the orthogonal to the hypersurfaces t = const, direction. 

The particular case with 77 = 4 corresponds to the system derived in [ pl]| . As we will show shortly, the latter does, 
however, not seem to admit a symmetric hyperbolic formulation for the evolution of the constraints. 



In pO[ the following symmetric hyperbolic system 1 20 
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B. The evolution of the constraints 



In order to solve Einstein's vacuum equations, one has to supplement the evolution equations (|4[|^,^) with the 
following constraints: 

0= C=y''''d''idabk-dkab)+l.O., 

0= Cj^d^Kaj-g'^'djKab + l.o., 

= Ckij = dkij — dkSij , 

= Cikij = d[idi;]ij ■ 

where C, C'i are the Hamiltonian and momentum constraints, respectively. 

Using the evolution system fl,|5|, |6|) fo r the main variables, one obtains the following principal part for the evolution 
system for the constraint variables ]10(|: 

doC = ^d'^Ck+l.o., (9) 

doa = ^—^d.,c-d''ci,,-d''CkJ + i.o., (10) 

V 

doCktj = I.O., (11) 

7] — A 

doCikij = ?7 9[,gi,](jCj) + gy-9[,Cfe] + Z.o. . (12) 

A system is symmetrizable if we can find a transformation that brings the principal part into symmetric form. In 
order to investigate under which conditions this system is symmetrizable hyperbolic, we split Cikij in its trace and 
trace-less parts: 

Cikij = Eikij + - {gi{iBj)k — 9k{iBj)i) + - gijWik , 

where Eikij is trace-less with respect to all pair of indices and where in terms of the traces Ski = C^^j.^^^, Aki = C^j/j^j^ 
and Vik = Cik/, 

4 12 4 9 12 

Bki — t; Ski — ^ Aki — -pVki , Wik — -Vik + — Aik ■ 
3 5 5 5 5 

Rewriting the principal part in terms of these variables, we find 

d^c - la'^Cfe + Lo., (13) 



9oQ = ^^d,C -^''Sk^~^''Ak^-^''Vk^+l.O., (14) 

■n 

^oSk^ = - Y (^a(feQ) - i gk^^'C,^ + l-o. , (15) 

^oAk^ = (1 - 77) Lo. , (16) 

doVk, = (^^-3^d[kC,]+l.o., (17) 

OqCmj = l.o.. (18) 
^oElk^J = l.o.. (19) 

Having split all the variables into their trace and trace-less parts the only natural transformations that remain are 
rescaling of the variables or linear combinations of the antisymmetric symbols Aki and Vki. First, looking at the terms 
depending on d in the evolution equation for C and vice-versa, we see that < 2 is needed in order to make the 
corresponding block in the principal part symmetric by a rescaling of C. Next, looking at the terms involving Ski in 
the evolution equation for Ci and vice- versa, through a similar reasoning, we obtain the condition > 0. Finally, if 
< r/ < 2, it is easy to see that the transformation 

Ak^ = Ak^ + Vm , Vk^ = (7?7/4 - 3)Afc, -f (?7 - 1)14, , 



brings the corresponding block into manifestly symmetric form. 
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We can summarize the result as follows: If < < 2, the principal part of the system ( |13| - |19| ) is symmetric with 
respect to the inner product associated with 



({/, U) = CC + C^a + ^ S'^'Sm + ^-V + V'^'Vm + C'^'Cmj + E^^^'EiMj , (20) 

1]^ 67] 8 — 67] 



where U = (C, Ci, Ski, ^m, Vki, Eikij)^ ■ It is interesting to notice that in a numerical empirical search for a value of 
T] that improves the stability of a single black hole evolution, Kidder, Scheel and Teukolsky found the value rj — 4/33, 
which lies inside the interval < ry < 2 On the other hand, the evolution of the original EC (77 = 4), for which 
we were not able to find a symmetrizer, according to does not perform very well in 3D black hole evolutions. Also 
of interest is that in the recent work of Lindblom and Scheel [T^ they note that the previously mentioned range of 
is also preferred. (See figure 5 of their paper; the range < 77 < 2 translates to —00 < 7 < —0.5). In that work they 
also study the dependence on another parameter z, which corresponds to a rescaling of a variable and therefore its 
effects do not influence the principal part of the equations, which is what determines the level of hyperbolicity of the 
system. 



IV. CHARACTERISTIC VARIABLES 

Here we discuss the characteristic variables of the main evolution system (HjsH) and of the evolution system ([l3[|ig| ) 
for the constraints. The characteristic variables are needed in order to give the boundary conditions of type (||), which 
yield a well-posed initial-boundary value problem. 

From now on we will concentrate on linearized (around Minkowski spacetime in Cartesian coordinates) gravity for 
simplicity. That is, the background metric is 

ds^ = —dt^ + 6ijdx^dx^ , 

where 5ij is the flat space metric. We also assume that our perturbations have vanishing shift and vanishing linearized 
dcnsitized lapse. In these cases, all lower order terms vanish (with the obvious exception of the right hand side (RHS) 
of (^)) and the evolution equations simplify considerably. 

We also choose our spatial domain to be a box Q. = [xmin,Xmax\ x [ymin,Vniax] x [zmin, Zjnax], though the analysis 
below could be generalized to other domains. 



A. Characteristic variables for the main system 



When linearized around flat spacetime, the main evolution system reduces to 

dtgi] = -'2.Kij , 

dtK,, = -^''fk^J , (21) 
dtfkij = -dkKij . (22) 

Since g.ij does not appear in the evolution equations for Kij and fkij , we do not consider its evolution equation in the 
following. Therefore, the system we consider has the simple form dtu = A^djU, with u — {Kij, fkij)'^ ■ 

The characteristic variables with respect to a direction n* are variables with respect to which the symbol P{n) = 
A^rij is diagonal. They can be obtained by first finding a complete set, ei, ...,624, of eigenvectors of the symbol (the 
ej's are called characteristic modes) and then expanding u with respect to these vectors. The coefficients in this 
expansion are called the characteristic variables. 

For the above system there are six characteristic variables with speed 1, six with speed —1 and twelve with zero 
speed. They are given by 



r- 



= Kij - fmj , (speed -1-1) 



'"ij ' = + /™i ' (speed -1) (23) 
^ti] = /fey ~ nk!mj , (speed 0) 

where = fkijn'^- In terms of these variables, the evolution system becomes 

dtvff' = ±d,^ - 5^'^^ , (24) 
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where 9„ denotes the derivative in the direction of n and d'^ = dk — n.fc9„ are the derivatives with respect to the 
directions that are orthogonal to n. 

1. Gauge, physical, and constraint-violating modes 

Before we proceed, it is interesting to give the following interpretation to the characteristic variables of the main 
system with respect to a fixed direction n*^: Consider a Fourier mode of u with wave- vector along , i.e. assume 
that the spatial dependence of u has the form e"'j^\ In this case, the linearized constraints assume the form 
L{n')u — 0, where L{n^) is a constant matrix. Also, since in the case we are considering the non principal terms 
vanish, characteristic modes of this form solve the evolution equations. 

Now we can check which characteristic modes (or combination thereof) satisfy the constraint equations and which 
do not. It turns out that all modes violate the constraints, except for the ones corresponding to the characteristic 
variables v'i^ and = v^^g — 6abv^(}^^ /2, where A, B denote directions which are orthogonal to n. 

Next, consider an infinitesimal coordinate transformation a;^ i-^ + X^. With respect to it, 

9ij ^ 9tj + 29(iXj) . 
Assuming that X* and Xj are proportional to e^"^^^ and using (^) we find that 

while the remaining variables are invariant. In particular, we can gauge away the variables Vm and fin''- The 
characteristic variables are gauge-invariant and satisfy the constraints. 
This suggests the following classification: We call 

• the variables gauge variables. 

• the variables v^^g physical variables. 

• and the remaining variables f^^^'' and v^i^ constraint violating variables. 

We stress that this classification is exact only if all the fields depend on the variable n^x'^ and the time coordinate t 
only. (In particular, it is exact in 1 -|- 1 dimensions since then there is only one space dimension.) Nevertheless, this 
classification sheds light on the boundary treatment below: For example, we will see that giving boundary data that 
is consistent with the constraints will fix a combination of the in- and outgoing constraint violating variables while 
we will be free to choose any data for some combination of the in- and outgoing gauge and physical variables (see 
equation ( p8| ) below). 

B. Characteristic variables for the constraints 

For our purpose, it is sufficient to find the constraint characteristic variables that have non-zero speed, since these 
are the ones that enter the boundary condition (^). In the system considered here there are three characteristic 
variables with speed 1 and three with speed —1, 

V^^^ =-a±^^^Cn,±{S,,,+A^,), (26) 

while the remaining variables have zero speed. 

V. CONSTRAINT-PRESERVING BOUNDARY CONDITIONS 



Here we start with boundary conditions for the constraints that ensure that they are preserved though evolution 
and then translate them into the boundary conditions for the main system. In order to do so we write the in- and 
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outgoing characteristic variables of the constraints in terms of (derivatives of) the characteristic variables of the main 



system. Well-posedness of the resulting constraint system is shown in the Section VI . 

In order to preserve the constraints we impose boundary conditions for their evolution of the form (H) with 5 = 
{g ^ would not preserve the constraints). That is, 

v}^^^L/vl-\ (27) 

where the matrix L is constant. As discussed in Sect. ||, the coupling matrix L must be "small enough" if we want 
to obtain a useful energy estimate. We will anaWze this question in the Section VI 

Next, our task is to translate the conditions ( |2?| ) into conditions on the main variables. Using the definition of the 
constraints, Eqs. (||), we find 

C = I , (28) 

Q = d^K.j - d,K , (29) 

Sm + Am = d'fmi ~ dj%,, + I [n^d'vi - 3d,vn] , (30) 



where Vk — f^.^^ — f^^k- Using Eqs. ( |28| , p9| , p0[ ) and the definition of the characteristic variables, Eqs. ( p3| , p6D , one 
obtains 



Vj^^ = d^.fl - a-.W ± i (l - ^ ) a- (24".^ - 2.a - + , (31) 



where the capital indices ^, _B = 1, 2 refer to transverse directions (i.e. directions orthogonal to n), and where we sum 
over equal indices. 

A first problem that arises in the above expressions is that with maximally dissipative boundary conditions one 
does not control normal derivatives 9„ of the incoming characteristic variables at the boundary and, therefore, cannot 
impose the above conditions, Eqs. (|l],^2|). However, since in Eqs. (|l|,^2|) the normal derivatives are present only on 
variables with speed ±1, we can use the equations (^4|) to trade normal derivatives by time and transverse derivatives. 
Doing so one obtains 



The second problem is to show well-posedness for the evolution system of the main variables. Conditions (|27^, 
together with Eqs. (^,^J) do not directly translate into conditions of the form (^) for the main variables since 
transverse derivatives appear in the expressions ( |3^ , |3^ ) pT| . In order to get a well-posed initial-boundary problem 
we look for appropriate linear combinations g ~ — ^ of in- and outgoing main variables and an appropriate 
coupling matrix L such that the conditions (^) with (^,^4|) can be incorporated in a closed set of evolution system 
at the boundary for the variables g and some variables with zero speed. This is discussed next. 

A. Neumann boundary conditions: evolution system on each face 

Consider the combinations dij — w-^^ — v^^^ and = v-^^ -I- u-^^^ (dy stands for difference |^, and Sy for sum). 
Noticing that 

yJl^^ + yy^ = dtdBB - d^SAn , 

+ V]"' = dtSAn + d^'dAB - dAdBB - 9^^™ + 29^1;^'^^ + j^i - Oa (2v^bL + dBs) , 
we see that one way of imposing boundary conditions such that the constraints are satisfied is through 

= v;(+) + kP, (35) 

= + Vj"^ . (36) 
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Equations amount to giving as boundary conditions for the constraints' propagation a coupling between the 

incoming and outgoing characteristic constraint variables as in Eq. ( p7| ) with (Lf) = diag{—l, 1, 1). These equations 
can also be seen as equations for dsB and SnA at the n-face, where dAB and dnn are a priori prescribed functions (i.e. 
they are not fixed by the constraints' treatment). Here dAB is defined as the traceless part of dAB- That is, 



J'AB 



^AB 



SABd(f 



where Sab = 1 if A = B and zero otherwise. 

In fact, Eqs. ( 35]3^ ) do not not contain as dynamical variables only dBB and SnA, but also some zero speed 

variables. Therefore, in order to get a closed system, we need evolution equations for the zero speed variables f^5„. 
These equations can be obtained from the evolution system of the main variables, equation (p5|): 



- a. JO) 



= dtv 



t"ABr. 



- BaSuB 



(37) 



Below we will show that equations (35 3qj37| ) constitute a symmetrizable hyperbolic system of evolution equations at 

the face orthogonal to n for the variables [dBB, SAn,v^ABn)^ where (Iab — dAB ~ 1/2 (Jab^cc and c?„„ can be freely 
prescribed. Since the domain we are interested in is a box, Eqs. (^^ 37) need to be evolved at each face. In order 
to do so one needs boundary conditions at each edge (intersection of two faces) . Below we will also explain how these 
boundary conditions are naturally fixed by compatibility conditions. 

Assuming one already has the solution to Eqs. (^,^,^^ at each face, the boundary conditions for the main 
variables, which are of the required form (y), are the following: 



„(+) _^(-) 

"BB ^ '^BB 



IBB 



,(+) - „(-) 



"nA 



SnA 



^nn nn 



where ds b and SnA are obtained from the evolution system 
variable dnn and the physical variable dAB freely. Note that 



dnn , = + dAB , (38) 

Tn at the n-face and where we can specify the gauge 



'2fn 



dAi 



-2fnAB + SABfnCC 



(39) 



In view of equation (^ and the fact that if the constraints are satisfied we have dkij — dkgij, we see that these are 
Neumann conditions on some components of the three-metric. 

Now we go back to the evolution system ( p5[]37| ) at the face and look at it in detail. In order to make the notation 

more compact we write d^"-* = dsB , s^"* = sati , h'^^B ~ ''^ABn ^'-'^ ^^'^ variables associated with the face orthogonal to 
n. Then we have 



^th^AB = -\dAS^B^ 



4-37? 



2-3?7 



dAdtJ, 



(40) 
(41) 

(42) 



where /i'"' denotes the trace of h^^g. One can check that as long as ?y ^ 2 this system is symmetric hyperbolic with 
respect to the inner product associated with 



4-3?7 



(43) 



where B — {d, sa, hAB) and Hab denotes the trace-less part of Hab, ^ab — hAB — {6ABh(f)/2. 



B. Neumann boundary conditions: Compatibility conditions at the edges of the faces 

The system we just introduced is defined on each of the six n-faces. The faces themselves have boundaries: the 
edges that join them. We need to ensure that the system on each face is well posed taking into account the boundary 
conditions at the edges. Since each edge is shared by two faces, this will translate into compatibility conditions among 
the various systems. 
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The system of interest has, at each edge, two ingoing modes with speeds vS/^ and 1 and two outgoing modes with 
speeds —^JiJ2 and —1. At face n, the corresponding characteristic variables with respect to a direction m^, are 



/372 2 V 3 " 3 



4-37? 



+ 



2-3?7 



w±i 



mp ' 



where p is a transverse direction orthogonal to m (and n). 

In terms of the original variables, the variables /i^r, dm},Si l on the n-face are 



BB 



2X„ 



MS 



— _2f 



Now let {n,m,p} be orthogonal directions naturally associated with the box (say, the standard Cartesian coordinates 
{x, y, z}). If we compare these variables in two different directions n and m, say, we obtain the following compatibility 
conditions: 



-2/;(") 



"'pn 



pm pn 

On the other hand, we have, by definition of the characteristic variables. 



(44) 
(45) 
(46) 



„(») 



+ ^3/2 -V3/2 V 3 ' " ' 

Therefore, the correct boundary conditions at the edges are 



in) 
+ 1 



W 



(") 



'4ft,(") 

^"'mp 



W 



(n) 



/372 

(n) 
''+1 



in) 
„(") 



2dl 



(47) 
(48) 



Notice that up to now the quantities that were freely specified were c?„„ and dAB one each of the n-faces. Therefore, 
in order to fix the boundary data for the systems that are defined on these faces we also have to a priori specify the 
quantities sL™^ — 2Knm — Sm"^ at the edges defined by the intersection of faces n and m. 



Imposing these boundary conditions automatically implies that the compatibility conditions 
satisfied. On the other hand, we have at the edge joining the faces n and m, 



and 



(h^n) _ ^(m)\ _ _\q ( in) _ (m)\ _ q 



where the last equality follows from imposing the boundary conditions at the edges. Therefore, these boundary 
conditions also imply that the compatibility condition ( ^6| ) is satisfied through evolution provided it does so initially. 



C. Dirichlet boundary conditions 

In a similar way to the Neumann case, one can obtain a closed system at the boundary for the variables 

/ J (0) (0) (0) \ 1 

(sBS,rfAn,'yVm.«ABS''f'sBA) requirmg 

= V^+^~VyK (49) 
0^vi+^+vi-\ (50) 

where 

-V^^ - = dtdAn + d^SAB - OaSBB - dASnn ■ 
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One also has to take into account the evolution equations for the zero speed variables: 

^^dtv% + \dAS,, . (51) 

In this case, one can freely specify s„„ and sab, which corresponds to Dirichlet conditions on some components of 
the extrinsic curvature: 

Snn = 2i4:„„ , SAB = '^K AB ~ ^AbKcC ■ (52) 

Assuming one already has the solution to Eqs. (49 j5^) at each face, the boundary conditions for the main variables, 
which are of the required form (|^), are the following: 

= -^is + SBB , v')^J = V^^J + dAn , vit^ = -vi~} + Snn , ^AB = "^AB + ^ AB ■ (53) 



^BB 



In terms of the variables s = sbb, dA — dAn, hA — Bry v^abb + (4 ^ '^v)i'>^BBA ~ ^Ann)' '^^^ rewrite the boundary 
system (|9[]5^) as 

dts = ^^d^dA-ld^hA, (54) 

dtdA = ^ dAS + dASnn " d^SBA , (55) 
dthA ^ -^^^^dAS + ^^^ (dASnn-d^SBA) ■ (56) 

This system is symmetrizable hyperbolic with respect to the inner product associated with 

A ^-iV,A\(u 4-377 



(B, B) = 4s^ + M'^dA + [h'' ^ rf^ j \ hA ^ dAj (57) 

where B — (s, dA, hA)- 

There is one ingoing and one outgoing mode, with speeds ±-y/3/2, respectively. The corresponding variables in a 
direction are 

r— , 1 /8- 377 , 

'^±7^ = V 3/2 s ± - (^^— d„ - 
We now turn to the compatibility at the edges of the faces. In terms of the original variables one has 

s(") = 2KbB , d^2^ = -2fnnA , h^A^ = 3?? /abb + (4 - 377) (/bBA - fAnn), 

= 2Knn , s^2l = '^Kab - SabKcc ■ 
From these expressions one can see that the only compatibility conditions are 

2 ^s^j^ — Stoot) , (58) 

^mm ~l" 2s^j„^ S^J + 2s[y2n ■ (59) 

Equation ( ^8|) fixes the boundary data for the ingoing variable: 

^W575 = + 4 V372(.(™) - 41)- (60) 



This condition will be used in the energy estimate (|6^) in order to prove well-posedness for the system at each face 
in the Dirichlet case. On the other hand, Eq. ( |59| ) is a compatibility condition at the intersection of faces n and m 
for the free boundary data. 



VI. WELL POSEDNESS 

We start by showing that the conditions (|35| , |36| ) and (^ , ^ ), for the Neumann and Dirichlet case, respectively, 
imply that the evolution system for the constraints is well-posed. Then we derive energy estimates for the closed 
system of evolution equations at each face, and using these estimates we show that the evolution of the main system 
is well-posed as well. 

Since we have already shown that all the evolution equations involved in our constraint-preserving treatment are 
symmetric hyperbolic, and since we have already cast all boundary conditions in maximally dissipative form, the 
main purpose of this section is to explicitly show that the different couplings are "small enough" with respect to the 
corresponding symmetrizers, in the sense discussed in Section O. 
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A. Constraint propagation 

Here we derive an estimate for the growth of the energy 

i=„„™.„..^/(t/.f/)A, 

Jn 

where {U, U) is defined in (pO|). Taking a time derivative and using equations (p^|T9|) we obtain, after integrations by 
parts, 



-2 

7, ^constraints — ^ 

at 



an 



4 - 2?7 



V 



da, 



where n*^ is the unit outward normal to the boundary dft of the domain. Expressing the integrand in terms of 
characteristic variables defined in (|2^), we obtain 



4 E..,„..,rr„.nts = \ i <5^M > V}^ > - Vr'V. 



where the last equation follows from the conditions (35,3g) and ( [49| , pOD in the Neumann and Dirichlet cases, respec- 
tively. Therefore, the initial-boundary value problem for the constraints is well-posed. In particular, this implies that 
zero initial data for the constraints implies that the constraints are zero at later times as well. 

B. Face systems 

1. Neumann conditions 

In order to show well-posedness for each system defined on the n-face r„ we consider the corresponding energy 
norm for the Neumann case, 



E(T^,N)= {B,B)da, 



where {B, B) is defined in (|4^). Taking a time derivative and using the above evolution equations we obtain 



d_ 



EiT„,N) = j^^ (wl^ - ^^-./iT^) + ^ ^ ""-i) + (^'^'^"" " d^'dAB) da. (61 



We use the boundary conditions ( p7p^ ), with si™'' = and dl™'' = for the moment, to get rid of the first term on 
the RHS of Eq. (p|). Using Schwarz's inequality, the second term is estimated as follows: 



4 / (dAdnn - d'^dAB) da < 4 f E^'^ 



(r„,w) 



where 



f^j^ {dAdnn - d^'dAs) (d^ d^n " dcd^"") da. 
Therefore, we end up with the estimate 

E^r^,N){tf'^ <E^r^,N){QY'^ + '2 I f{s)ds (62) 



and see that the energy is bounded and the bound is determined by the norm / of the free data on the boundary. 

The assumptions si™'' = = dirp' can be easily relaxed with an argument similar to the one we presented in the 
paragraph following equation (^. That is, one defines a new variable that incorporates the inhomogeneity in the 
energy generated by the term that appears in the case of a non-smooth boundary and obtains an energy estimate for 
the new variable, since the variable redefinition is finite and well defined. Well posedness follows immediately. 
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2. Dirichlet conditions 
The Dirichlet case is similar to the Neuniann one. The energy is now given by 

E(r„,D)^ {B,B)da, 
with {B,B) given by Eq. (|57|). We obtain the estimate 



(63) 



where 



J {dASnn - d^SAs) {d^Snn " dcs^^) da, 



and we can proceed as in the Neumann case, using in this case the boundary condition (|60|). 

C. Main system 

Having obtained a bound for each closed system defined on a face we can obtain a bound for the main evolution 
variables by the standard techniques described in section II where we now have the boundary conditions ( |38| ) and 
( ^ ) for the Neumann and Dirichlet cases, respectively, which are of the required form (^). 

VII. SUMMARY 

The well posed, constraint preserving boundary conditions presented in this paper can be summarized as follows. 

A. Neumann case 

• Free data: 

At each face (say, the n-one), the three quantities d„„ and dAB must be a priori defined (subject to standard 
compatibility conditions with the initial data). These quantities correspond to the normal (that is, in the n 
direction) derivative of the normal and the transverse, traceless parts of the three metric, see Eq. (p9|). One of 



these variables is gauge and the other two are physical, in the sense discussed in Section IV. In addition, the 



quantities Sm'' = 2iir„,„ should be prescribed at each edge defined by the intersection of faces n and m. These 
quantities must satisfy the compatibility conditions s^™^ = ^ . 

• Evolution systems on faces: 

The 2D symmetric hyperbolic 7x7 system (^ , ^ , ^2| ) is evolved on each face. This system needs boundary 
conditions at the edges of the corresponding face. They are given by equations (^, ^). 

The solution to each of these systems provides the three quantities d'^^ and s^"^ at each of the six n-faces. 

• Main evolution system: 

The main system (21 2^) is evolved in the 31? domain. This system needs boundary conditions, at each face. 



for the six incoming characteristic modes. These boundary conditions at, say face n, are given by Eq. (|3^ 
where the needed information for three of these boundary conditions is provided by the a priori specified d; 
and dAB, while the other three are given by and s^^J^. 
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B. Dirichlet case 

Free data: 

At each face (say, the n one) , now the three quantities s„„ and sab must be a priori given. They correspond to 
the time derivative of the normal and transversal, traceless part of the three metric, see Eq. (p2|). These three 
quantities have to satisfy the standard compatibility conditions with the initial data, but also some compatibility 
conditions at edges, Eqs. (^,^). They are also gauge and physical variables, in the sense discussed in Section 



[V 



• Evolution systems on faces: 

The 2D symmetric hyperbolic 5x5 system ( |5^,|55| , p6| ) is evolved on each face. This system needs boundary 
conditions at each edge. They are given by Eq. (pOf). 

The solution to each of these systems provides the three quantities s^^ and at each of the six n-faces. 

• Main evolution system: 

The main system (21 2^) is evolved in the 3-D domain. This system needs boundary conditions, at each face. 



for the six incoming characteristic modes. These boundary conditions are given by (^3|), where the needed 
information in three of these boundary conditions is provided by the a priori specified s„„ and sab , while the 
other three are given by s^'^ and c^An- 

VIII. CONCLUSIONS: 

We have studied the system of evolution equations for the constraints for a subfamily of the generalized Einstein- 
Christoffel symmetric hyperbolic system. We have shown how to give boundary data for the constraints in such a way 
that it translates into boundary data for the main system that yields a well posed problem, both for the main system 
and the system of evolution equations for the constraints. We have studied the case of a boundary that is not smooth, 
as is the case of the usual cubic boxes used in numerical relativity. This required additional care at the boundaries 
of each face, with the ensuing compatibility conditions. It should be noted that the energy estimates derived do 
not necessarily guarantee the existence of a smooth solution in the presence of corners even with the compatibility 
conditions we presented. Further work is needed to establish smoothness of the solution. 

Our analysis was carried out for the case of linearized gravity around Minkowski space-time. It is expected that 
similar techniques will be useful in the case of other background space-times and also in the non-linear case. We 
will discuss these generalizations in future papers. Also, since we have followed a systematic approach and have not 
taken any advantage of the gauge choice, in principle it should be possible to apply the same procedure to symmetric 
hyperbolic formulations with live gauges [l^ . 

We have also found that, at least with the formulation of Einstein's equations here used (the generalized EC), 
the Neumann and Dirichlet cases are in fact the only ones for which well posedness can be established through the 
techniques used in this paper (see the appendix). More specifically, we have found that these two cases are the only 
ones in which closed systems at the faces can be obtained. However, this does not mean that these are the only well 
posed cases, since in the initial-boundary value problem an energy estimate is a sufficient but not necessary condition 
for well posedness. Also, giving the appropriate ((in„, ^("^ab, Snm) data in the Neumann case, or the appropriate 
{srm,SAB) in the Dirichlet case, one should be able to recover any spacetime in any slicing. This is because our 
constraint-preserving treatment makes sure that one is solving not only Einstein's evolution equations but also the 
constraints. But it does not make any restriction on the space of solutions to the Einstein equations. However, it is 
not clear how to choose these "appropriate" boundary conditions, without any a priori knowledge of the solutions, in 
order to model an isolated source, given that the boundaries are at a finite distance. This same problem appears in 
similar approaches |l . One possible solution is to provide these functions through Cauchy-characteristic M, |l^ or 



Cauchy-perturbative matching |18|, or to resort arguments using the peeling property. Another possibility would be 
to impose a "no incoming radiation" condition. However, it is not clear how to do this within formulations that have 
as extra variables first but not second spatial derivatives of the three metric. This might be remedied by repeating 
this construction for other systems of evolution equations where the variables that represent gravitational radiation at 
the boundary play a more central role, as in For instance, a system of evolution equations of higher order where 
the Weyl tensor is the fundamental variable, would be suitable for this purpose. 

With the results of this paper one can now assure that both the evolution equation for the main variables of the 
problem and the evolution equations for the constraints are well posed on a manifold with (non smooth) boundary. 
This allows to evolve initial data that satisfy the constraints beyond their domain of dependence, as is of interest in 
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numerical simulations of the binary black hole problem. Moreover, well posedness opens the possibility of constructing 
numerical schemes for which numerical stability can be rigorously proved. 
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APPENDIX A: CLOSING THE BOUNDARY SYSTEM 

In Section we showed how to construct a closed system at the boundary by taking appropriate linear combinations 
of characteristic variables. We also chose a particular combination of ingoing and outgoing constraints and showed 
that it gives rise to a system of partial differential equations that lives on the boundary. To close this system we had 
to include the evolution of some zero speed modes. A question that arises is whether there other ways of closing the 
boundary system, apart from the Neumann ( ^o[]4^ ) and the Dirichlet ( p^ - |5^ ) case. 

To answer this question we will make no assumptions on the coupling matrices R and L. The boundary condition 
for the system of the constraints is assumed to be 

V}+^ = L^vj'^ , (Al) 

where Li^ is a 3 x 3 coupling matrix and V^^'' is given in (|3^ ) and (0). At the boundary data must be given to the 
ingoing modes. We will assume that they satisfy 

vlt^ =Rv'''vi;^+b.,, (A2) 



where is a 6 x 6 coupling matrix and bij is the boundary data. If we insert (A2) into (Al) we obtain a system 



which contains derivatives of the boundary data bij , of the outgoing modes v^j ^ , and of the zero speed modes v!^^j . 
This system can be solved for the time derivatives of three of the bij, namely bsB and bnA- To the remaing bij one 
can give arbitrary data and consider them as source terms. In order to close the system the coefficients that multiply 
terms containing outgoing modes ^ must vanish. After imposing this condition the coupling matrices R and L, 
which in general depend on 45 parameters, depend on one parameter only (apart from rj). The zero speed modes that 
appear in the rhs of this system cannot be eliminated by any choice of this parameter. Therefore in order to close the 
system one has to enlarge it by including the evolution of at least the zero speed modes that appear in the rhs. The 
requirement that the evolution of these zero speed modes do not contai n an y spati al d erivatives of outgoing modes. 



forces the couplings R and L to be the ones the we used in subsections ( |V A| ) and ( |V C| ). 

Summarizing, the Neumann and the Dirichlet cases are the only ways that one can obtain a closed system at 
the boundary. Furthermore, as we have shown in this paper, the boundary system is symmetric hyperbolic and the 
coupling matrices R and L are "not too large". Any other choice of coupling matrices would lead to a system for 
which the techniques used in this paper to prove well-posedness cannot be applied. 
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This system corresponds to System 3 of [nof with their parameter z set to zero. 

Note that this difficulty does not arise in 1 + 1 dimensions since then there are no transverse derivatives. In that case one 



can simply set the ingoing constraint violating variables «^ and ii^^ 



„(+) 



to zero. 



Notice that d is used for this quantity as well as for dkij 
since both objects have different number of indexes. 



dkQij, equation ([^). There does not seem to be risk of confusion 



